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Abstract 

In this work we study the dynamical properties of a finite array of nanomagnets in artificial 
kagome spin ice at room temperature. The dynamic response of the array of nanomagnets is studied 
by implementing a "frustrated celular automata" (FCA), based in the charge model and dipolar 
model. The FCA simulations, allow us to study in real-time and deterministic way, the dynamic of 
the system, with minimal computational resource. The update function is defined according to the 
coordination number of vertices in the system. Our results show that for a set geometric parameters 
of the array of nanomagnets, the system exhibits high density of Dirac strings and high density 
emergent magnetic monopoles. A study of the effect of disorder in the arrangement of nanomagnets 
is incorporated in this work. 
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I. INTRODUCTION 



Based on Dirac's prediction about the existence of magnetic monopoles, considering the 
quantization of the electrical charge [1], many research teams have attempted to find these 
monopoles [2], but their results have been inconclusive. Attention is currently directed at the 
study of elemental excitations in natural and artificial spin ice systems [3-6]. Excitations 
appear in natural systems that are equivalent to the magnetic monopoles connected by 
Dirac strings, which can be observed indirectly [1, 7-9]. In artificial systems a pattern of 
magnetic nanoislands is created with lithographic techniques (with a square lattice or kagome 
lattice) equivalent to natural spin ice systems [10]. In the case of artificial spin ice systems, 
the nanoisland are arranged such that there is frustration to minimize dipolar interaction. 
This allows for designing geometries where emergent excitations are present, equivalent 
to the magnetic monopoles of natural spin ice systems. As well, with the artificial spin 
ice systems it is possible to visualize the configuration of magnetic moments directly with 
different microscopic techniques. This allows for their study at ambient temperature and the 
simulation of the behavior of natural spin ice systems [10-16]. Mengotti et al [17] recently 
published the results of direct observations of emergent monopoles and associated Dirac 
strings. These results are the first direct confirmation of the reduction in the dimensionality 
of the system as a result of the frustration in the artificial spin ice system. Mengotti's work 
provides experimental evidence of reverse magnetization through the Dirac string avalanches 
that join monopole-antimonopole pairs. The same work provides results of Monte Carlo 
simulations that are in agreement with experimental results. The excellent work of Mengotti 
et al [17] provides very valuable information about the dynamic of these elemental excitations 
and about nucleation and avalanches in reverse magnetization. Equally, many questions 
were raised that have lead to studying some aspects of the studied system by the authors 
of this work [17]. For example, it is expected that the formation of monopole-antimonopole 
pairs begins at the ends of the sample and that the strings advance toward the center. 
In the aforementioned report [17], only the central region of the array is shown, without 
the possibility of verifying the initial formation of the strings. Equally, the Monte Carlo 
simulations of these and other authors, considering infinite systems, are not very efficient 
when studying at the beginning of the reversion of magnetization. One the central points 
of this paper is to introduce the use of cellular automata to the field of frustrated magnets. 
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In this work, we study the effect of the aspect ratio and size of the system on the formation 
of Dirac strings. A study of the effect of disorder in the arrangement of nanomagnets is 
incorporated in this work. We considered the system as a finite array of magnetic nanoislands 
and took into account the possible disorder in the value of the nanoisland moments. The 
simulation is carried out using a determinist type of cellular automaton [18], which considers 
the coordination number of the elements of the grid. The interaction among the nanoislands 
is studied under the magnetic charge model and dipolar model. 

II. MAGNETIC NANOISLAND ARRAYS AND EMERGENT MONOPOLES 

The system studied in this work is a magnetic nanoisland array in a hexagonal lattice. 
Figure 1 shows a scheme of the magnetic bars arranged on the sides of a hexagon. Three 
nanomagnets converge in the vertices of the hexagons, as can be appreciated in Figure 1. 
Geometrically, we can define 2 non-equivalent vertices in the hexagonal lattice, vertices A and 
B, respectively. The two non-equivalent vertices form the unit cell of the entire array. The 
cell is shown in red in Figure 1. In the magnetic charge model, the charges are concentrated 
in these vertices. We define the magnetic charge in a vertex as -1, when two poles south 
and a north pole of the three nanomagnets that form the vertex converge. Equally, we 
define the magnetic charge as +1, when two poles north and a south pole converge in the 
considered vertex. Figure 2 shows the array with an applied magnetic field. The upper 
left of Figure 2 shows the array submitted to a field in the direction of the negative axis x. 
In the lower left of the Figure the field is shown directed to the right. The nanomagnets 
with the x component of the magnetic moment directed to the left are represented by gray 
squares, and the nanomagnets with the x component from the magnetic moment directed 
to the right are represented by dark blue squares. When the sample is totally magnetized 
in directions +x or —x, all the type A vertices have +1 (-1) charges, while all the B vertices 
have -1 (+1) charges. This can be seen in Figure 2. 

We suppose that we are under the condition of total magnetization, with the magnetic 
field directed to the left (upper left of Figure 2). Under this condition, we define a positive 
and mobile monopole, if a nanomagnet converging in a class A vertex inverts its magnetic 
moment. The charge of vertex A goes from qA = — 1 — > q\ — +1 =^ Ag^ = +2. If this is 
produced in a B type vertex, we define a negative monopole and would have qs = +1 — > q* B = 
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Figure 1: Scheme of the nanomagnet array 
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Figure 2: Scheme of the nanomagnet array when the sample is totally saturated by the influence 

of an external magnetic field. 

— 1 ==>- Aqs = —2, where q* A and q* B represent the charge of the vertex after the inversion 
and qA and q B represent the charges of the vertices A and B, respectively, in the their 
initial states. In this manner, when a nanomagnet inverts its magnetic moment, emerges 
a monopole-antimonopole pair. Figure 3a - 3d, shows the generation of a pair of mobile 
monopoles and their separation, giving rise to a Dirac string. If the three nanomagnets that 
converge in a vertex invert their moments, the condition Ag^ = +2 and Ag# = —2, is also 
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generated, but in this case the monopoles remain trapped and do not move through the 
sample (Figure 3e). 
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Figure 3: Scheme of the creation of a monopole-antimonopole pair, the associated Dirac string and 
trapped monopole. 



III. CELLULAR AUTOMATON, THE MAGNETIC CHARGE MODEL AND 
DIPOLAR MODEL. 

A cellular automaton (CA) is a mathematical structure used to model the dynamics of 
complex systems. It is formed by many simple entities that interact locally. A variety 
of models based on CA have been used to efficiently study problems in biology, physics, 
chemistry, engineering and material sciences [19-20]. They represent an excellent alterna- 
tive to models based on differential equations and to Monte Carlo algorithms because they 
can simulate highly complex systems with a low computational cost. The first attempt to 
use CA in the study of magnetism was the model proposed by Vichniac [21], which was 
subsequently developed by Pomeau [22] and Hermann [23] and is termed the VPH model. 
This is being used to resolve an Ising type spin system. To avoid a "feedback catastrophe", 
the automaton is updated in more than one step. The model functioned well at high tem- 
peratures (T > Tc), but failed at low temperatures. Subsequently, Ottavi et al [24] used a 
microcanonical algorithm in a CA to resolve the Ising spin system. A determinist version of 
this model provided acceptable results at a low temperature [24]. Owing to the popularity 
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of the different types of Monte Carlo algorithms used in problems associated with spins, 
the development of CA models for these systems did not continue. In this work we used a 
CA model, different from previous models, developed specifically to resolve the dynamic of 
frustrated spins in artificial spin ice systems. This model allows the simulation of spin ice 
systems efficiently [18]. 

A. Frustrated Cellular Automaton (FCA) 

This model was conceived for frustrated systems whose dynamic develops at zero tem- 
perature or the equivalent. In the case of artificial ice spin systems, each nanomagnet has a 
shape anisotropy with energy on the order of 10 4 K. This means that the thermal fluctua- 
tions in the configuration of moments are negligible at room temperature and consequently 
the system behaves like a system at zero Kelvin. To study this system in particular, we 
define the cells of the automaton in the vertices of the hexagonal structure. The ends of 
three nanomagnets converge in each vertex. Figure 4a shows the structure of the FCA. The 
automaton is updated as follows: 

1. A class of diagonal nanomagnets (dl or d2) is randomly selected at each stage of the 
algorithm. 

2. It covers all the vertices of the automaton. 

3. The moment of the nanomagnet of the chosen class is inverted and if the total energy 
decreases (and this difference is greater than the energy barrier for investment), the change 
is accepted. The energy term is explained in the next paragraph. This new configuration is 
maintained in an auxiliary array. 

4. The same is then done with the horizontal nanomagnet (H), not considering the 
nanomagnet of the non-selected class. 

6. The auxiliary configuration is copied in the definit configuration and we return to step 

1. 

In this way, we study the dynamic of emergent monopoles with a deterministic model 
that allows studying finite systems, considering the effect of the edge and the size of the 
system on the dynamic. As well, it allows us to incorporate impurities in the sample and 
study the effect of these impurities in the monopoles and associated Dirac strings. 
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a) 




cells forming the FCA 



b) 



charge model 



Figure 4: a) Scheme of the FCA cells, b) Scheme of replacing the magnetic moments by the 

magnetic charge model. 

B. Magnetic charge model. 

The moment rrt of each nanoisland in this model is replaced by two charges (one positive 
and the other negative), located at the ends of the nanomagnet, as shown in Figure 4b. 
The magnitude of each charge is q = y, where / is the length of the bar. The total charge 
in each vertex is the sum of the three charges associated with the vertex. Vertex j gives 
Qj = '^^Qk- The total energy of the system is given by the expression: 



2 4tt „ ' L l J 



V { y 113 (1) 

fi, i = j 

The term for i ^ j takes into account the interaction among the vertices of the array. The 
term i = j considers the energy of the site. This term considers the interaction among the 
ends of the three nanomagnets that converge. Figure 5 shows a scheme of the configuration 
of charges in each vertex and the parameters of associated length. 





Figure 5: Scheme of the configuration of charges in each vertex and of the length parameters 
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The energy in each vertex is given by the expression: 



h ~4n\ d + d + d i {) 
In accordance with the parameters of the hexagonal lattice d = ^ (a — I). Defining 

2 

Qo = m h an d writing the energy in units of the total energy can be written as: 

[/ = l ij lJ (3) 

[ ^ {?i?2 + gigs + g2gs} , i = j 

With £ = 1 — l /a. When the automaton is updated, the change in total energy is registered 
(using equation 3), that is, under the magnetic charge model. The interaction of the charges 
with the applied magnetic field is added to this term and the anisotropy energy term . 



C. Dipolar interaction model 

The dipolar interaction between sites on the kagome lattice is given by 

U dip = ~^ Yl ^ { 3 ' ^) (™j -rn) -m-rhj} (4) 

where m ; = ajmej is the magnetic moment at site i. Here the pseudospin ctj = ±1 denotes 
the projection of the spin onto the anisotropy (local ising) directions e>j (directed along 
the links of the honeycomb lattice) and m is the magnetic moment of a nanoisland. The 
interaction of the charges with the applied magnetic field is added to this term and the 
anisotropy energy term . 



IV. SIMULATION AND THE RESULTS OBTAINED 



A. Monopoles and Dirac Strings with the charge model and dipolar model. 

The first system studied, is a sample of (50 /im x 50 /im) and 5,800 nanomagnets, with 
impurities. The lattice constant (the distance between two adjacent vertices) has a value of 
a = 577nm. The random magnetic moment of individual islands is given by m = mo/3, where 
(3 is a dimensionless Gaussian random variable with ((3) = 1 and s = (((/?— (P) 2 ) )) l ^ ■ I n 
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this first simulation, s = 0.13 and the energy barrier for the magnetic reversal is 98 x 10 4 kT, 
with T, the room temperature. This energy corresponds to a magnetic field of /j, H = 12, 3 
mT used in the simulation montecarlo [17]. The quantities N T and cr M shown in all the 
studies are represented in units of mo ^°' 6S and the magnetic field is shown as normalized to 
the coercive field He- First, we do the simulation using the dipolar model. Figure 6a shows 
the hysteresis curve, experimental data of the hysteresis curve [17] and the total density 
of monopoles. Figure 6b shows the density of mobile monopoles and experimental data in 
function of the magnetic field. The agreement with experiment is very good. The discrepancy 
in the size of the plateau is because in this region, the simulation is very sensitive to the 
sample. This sensitivity is also reflected in the Montecarlo simulation [17]. It is important 
to note that at the left end the sample ends with type "A" vertices and in the right end 
with type "B" vertices. The monopoles generated in the latter vertices, when the moments 
of the nanoislands are inverted are not mobile, and consequently cannot migrant through 
the sample. The latter implies that the pair generated in the ends only contributes a mobile 
monopole that is shifted toward the center of the sample because of the applied magnetic 
field. Figure 7 shows a scheme of the complete sample (all the magnets) and the region that 
is considered for the statistical analysis (red rectangle), for values of the magnetic field near 
the coercivity. This figure does not include the color of each vertex to better appreciate the 
mobile monopoles. The impurities present in the sample (which were deposited randomly) 
produce mobile monopoles pairs in every region of the surface and both the positive and 
negative pole can move. This is clearly illustrated in Figure 7, which simulates the dynamic 
of system. We can note that some pairs of emergent monopoles appears in the center of 
the sample. These emergent monopoles, which come from impurities and are generated in 
the central part of the nanomagnet arrays, move toward the ends, thus extending the Dirac 
string. 
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Figure 6. Dipolar model (fioH c = 18.05 mT). a) Hysteresis curve (experiment [17] and simulation) 
and total density of monopoles. b) Density of mobile monopoles (experiment [17] and simulation). 
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Figure 7. Scheme of the simulation of the dynamic of monopoles and Dirac strings for a system 
with 5,800 nanoislands with impurities (dipolar model). 



We now examine the behavior of the system using the magnetic charge model. All param- 
eters are identical to the dipolar case, but now we add the size of the magnets (required in 
this model). We use the value I = 430 nm, which allows us to obtain, the experimental en- 
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ergy of site. Figure 8a shows the hysteresis curve, experimental data of the hysteresis curve 
[17] and the total density of monopoles. Figure 8b shows the density of mobile monopoles 
and experimental data in function of the magnetic field. We can appreciate some impor- 
tant differences compared to dipolar simulation. First, the maximum density of magnetic 



monopoles reaches a value of 0.082 



monopoles 
sites ' 



compared to the value of 0.095 



monopoles 
sites 



111 



dipolar simulation. 

Moreover, the total density of monopoles presents a curve with a slope greater than in 
the case dipolar. The hysteresis curve has no plateau characteristic of experimental data 
and dipolar simulation. These differences in results are due to the energy of the site (charge 
model), which prevents the initial reversal of nanomagnets H (Figure 4). This behavior can 
be seen by comparing Figures 7 and 9. 

We can appreciate a good agreement between the statistical results of the dipolar simula- 
tion with experimental statistical results. However, the images (XMCD) of the experimental 
paper in the first phase of magnetic reversal, are more similar to the images obtained with 
the simulation of magnetic charges. 
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Figure 8. Magnetic charge model {[iqH c = 11.01 mT). a) Hysteresis curve (experiment [17] and 
simulation) and total density of monopoles. b) Density of mobile monopoles (experiment [17] and 
simulation) . 
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Figure 9. Scheme of the simulation of the dynamic of monopoles and Dirac strings for a system 
with 5,800 nanoislands with impurities (charge model). 

B. The effect of the size and aspect ration of the system on the dynamic of emer- 
gent magnetic monopoles 



We summarize our study of the density of mobile monopoles in function of the size of the 
system. For each size, the simulation was repeated 100 times. The error bars in the graphs 
correspond to standard error. Figure 10 shows the maximum density of mobile monopoles, 
in units of m °"°^°' es , using the dipolar model, for two concentrations of impurities s = 0.13 
and s = 0.25. All systems, shown in Figure 10, have square geometry. Axis x in Figure 10 
shows the side of the square in fim. From the Figure, we can appreciate that for all systems, 
the density of the mobile monopoles is higher for the case with S = 0.25. As the size of 
the system decreases, the density with impurities approaches to the maximum value of the 
m0 ™ite° leS ■ Our study verifies that the maximum density of magnetic monopoles, decreases 
exponentially with system size. Figure 11 shows the same studied, using charge model. We 
can verify a qualitative behavior similar to the dipole model, but with lower density values, 
(discrepancy explained above). 
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Figura 10. Dependence of the density of the monopoles on the size of the system. (Dipolar model). 
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Figure 11. Dependence of the density of the monopoles on the size of the system. (Charge model). 

Finally, we show the results of the study considering the aspect ratio of the system. Figure 
12 shows the maximum value of the density of mobile monopoles in function of the ratio j^, 
where Lx is the length of the sample in the direction parallel to the magnetic field and Ly 
is the direction perpendicular to the magnetic field. All systems studied have an impurity 
concentration, equivalent to s — 0.13. Our results show that the density of monopoles is 
highly sensitive to the aspect ratio of the system. 
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Figure 12. Dependence of the density of the monopoles on the aspect ratio of the system. 



V. CONCLUSIONS 



In this work we have studied the dynamic of magnetic monopoles emergent in an artificial 
spin ice system. An algorithm based on a frustrated cellular automaton was used in the 
magnetic charge model and dipolar model. The great advantage of the model is that it can 
make efficient simulations of highly complex phenomena in real time with a minimum of 
computational requirements. The model represent a perfect complement to the methods 
based on Montecarlo algorithms, to study the elemental physics of problems with classical 
and quantum entities. Our results show that the number of emerging magnetic monopoles, 
depend of the sample size, the aspect ratio and the concentration of impurities. This allows 
us to consider a possible engineering in the creation of these systems for applications, for 
example in information technologies [25]. 
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